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It is a well known fact that subdiffusion equations in terms of fractional derivatives can be ob- 
tained from Continuous Time Random Walk (CTRW) models with long-tailed waiting time distri- 
butions. Over the last years various authors have shown that extensions of such CTRW models 
incorporating reactive processes to the mesoscopic transport equations may lead to non-intuitive 
reaction-subdiffusion equations. In particular, one such equation has been recently derived for a 
subdiffusive random walker subject to a linear (first-order) death process. We take this equation as 
(— i ' a starting point to study the developmental biology key problem of morphogen gradient formation, 

both for the uniform case where the morphogen degradation rate coefficient (reactivity) is constant 
and for the non-uniform case (position-dependent reactivity). In the uniform case we obtain expo- 
nentially decreasing stationary concentration profiles and we study their robustness with respect to 
perturbations in the incoming morphogen flux. In the non-uniform case we find a rich phenomenol- 
ogy at the level of the stationary profiles. We conclude that the analytic form of the long-time 
morphogen concentration profiles is very sensitive to the spatial dependence of the reactivity and 
the specific value of the anomalous diffusion coefficient. 
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•<-J ' I. INTRODUCTION 

>> 

tIh ' Fractional diffusion equations are a powerful tool to study anomalous transport processes, i.e., processes in which 
the mean square displacement (x 2 ) of a randomly moving particle displays the long time-behavior (x 2 ) ~ K~t^ , where 
7 is the anomalous diffusion exponent and if 7 is the so-called anomalous diffusion coefficient. When < 7 < 1, one 
has sublincar growth of (x 2 ) (subdiffusion), while for 7 > 1 one speaks of superdiffusion. As it is well known, the 
^ ' classical diffusion equation corresponding to the 7=1 case can be obtained from an average over the trajectories of 
a Markovian random walk in the limit of large time scales and long displacements. In contrast, stochastic transport 
' processes governed by anomalous diffusion equations reflect memory effects at the microscopic level. In particular, 
one can show that a suitably defined non-Markovian hopping process, namely a Continuous Time Random Walk 
(CTRW) with a long-tailed waiting time distribution, yields a subdiffusion equation in terms of the Riemann-Liouville 
fractional derivative [lj . This fractional subdiffusion equation can be taken as a starting point to deal with a number 
of biologically relevant problems, e.g. the localization of a target protein by a sea of subdiffusively moving ligands in 
the intracellular environment [U, Q ■ In this case, the complexity of the cell medium results in the ligands encountering 
a large number of obstacles, barriers, etc. in the course of their trajectories. In the framework of CTRW models, the 
effect of this crowded environment can be partly captured using waiting time distributions; subsequent averaging of 
the resulting equations over trajectories then leads to the associated fractional subdiffusion equations. 

While anomalous diffusion and in particular subdiffusive processes play a central role in Nature as a manifestation of 
underlying memory effects at a microscopic level, the situation where the particles simultaneously undergo anomalous 
transport and reaction (understood as a particle creation, destruction or transformation process) is also very common 
and important from the point of view of biological applications. In the example of the target protein and the ligands 
given above, one could allow e.g. for the possibility of the ligands undergoing a degradation process as they sweep 
the cell medium. Degradation implies a change in chemical structure which results in the ligands losing their ability 
to interact with the target; for practical purposes this kind of transformation can therefore be regarded as a "death" 
or "evanescence" process. 

In what follows, we shall focus on yet another biological process where degradation/death plays a central role, 
namely morphogen gradient formation. The location, differentiation and fate of many embryonic cells is governed by 
the spatial distribution of special signaling molecules called morphogens. Standard models of morphogen gradient 
formation assume that a specific part of the embryo secrets morphogens at a constant rate. The secreted morphogens 
then undergo degradation as they disseminate through the tissue and a concentration gradient builds up. Different 
target genes in the embryonic cells are activated above different morphogen concentration thresholds, implying that 
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the cell response to the local environment will depend on how large the concentration is. Thanks to this differential 
response, cells are able to interpret the morphogen gradient and translate it into specific "code" for their further 
development via the expression of the relevant genes. 

Traditional models of morphogen gradient formation are based on classical diffusion equations with a linear degrada- 
tion term. Here, we aim to go one step further and allow for the possibility of anomalous transport, as memory effects 
are likely to strongly influence the stochastic motion of morphogens in the complex embryonic environment. A re- 
markable property of morphogen gradients is their robustness against changes or fluctuations in the rate of morphogen 
production or degradation. An interesting question which we aim to study is the interplay between subdiffusion and 
robustness with respect to such perturbations. 

If one accepts the idea that morphogens perform subdiffusive motion as a result of strong dispersion in their waiting 
times between consecutive jumps, great caution must be exercised when incorporating the effect of a simultaneous 
degradation process to the transport equations because of the non-Mar kovian character of the latter. Several recent 
works indeed illustrate that heuristic equations where one has separate terms for the reaction and the transport 
process may lead to unphysical results, e.g. negative particle concentrations (see e.g. @). Therefore, the derivation 
of physically correct (but not necessarily intuitive) reaction-subdiffusion equations calls for the use of an extended 
CTRW formalism where the effect of reaction is incorporated at a mesoscopic level of description. In a recent work 
the authors have shown by means of Fourier-Laplace techniques that CTRW models extended in such a way yield 
equations which (in addition to a standard, purely reactive term) display a mixed reaction-transport term containing 
both the reaction rate coefficient (reactivity) and a Riemann-Liouville fractional derivative with respect to time. 

The remainder of the paper is organized as follows. We first give a brief reminder of classical reaction diffusion 
equations used for modeling of morphogen gradients and subsequently discuss how to extend such equations to account 
for anomalous transport via fractional derivatives. We subsequently focus on the specific case of uniform reactivity and 
assess the robustness of the resulting stationary concentration profiles with respect to perturbations of the incoming 
morphogen flux. Next, we turn to the non-uniform case and discuss the long-time behaviour of the profiles for several 
specific situations, namely the case of a piecewise constant reactivity (which not always leads to a stationary profile) 
and the case of a decaying reactivity respectively given by an exponential and a power law. In some cases, we also 
provide numerical simulation results based on a CTRW model for evanescent particles and find excellent agreement 
with the analytic results obtained from the fractional diffusion equation approach. Finally, a summary of results and 
possible avenues for future research in this area are given in the conclusions section. 



II. CLASSICAL REACTION-DIFFUSION EQUATION WITH LINEAR DEGRADATION 

The cornerstone of many studies concerning morphogen gradients is the classical one-dimensional reaction-diffusion 
equation 



dc(x, t) 

at 



K 



d 2 c(x,t) 
dx 2 



k(x, t) c(x, t), 



(1) 



where the evolution of the concentration c(x, t) is described by a Fickian term (characterized by a classical diffusion 
coefficient K{) and a linear degradation term (characterized by the reactivity k(x, t)). Eq. ((T|) is then solved subject 
to the radiation-type boundary condition 



dc(x, t) 



dx 



dc(x, t) 



x=0+ 



dx 



■ Jo- 



(2) 



x=0- 



This boundary condition simply states that a constant flux of morphogens jo is injected at the origin x — 0. The 
simplest case is given by a constant degradation rate k(x,t) — k, which yields the exponentially decaying stationary 
profiles: 



c(x, oo) = J ° e~ V * N . (3) 

Despite its simplicity, the exponential dependence of Eq. © captures surprisingly well the rapid concentration 
decay displayed by real profiles. However, the separate determination of K\ and k poses significant experimental 
difficulties and in most cases only the characteristic length yTKyjk can be unambiguously determined. This opens 
the door to the possibility of considering non-Markovian generalizations of Eq. (JTJ) which are also compatible with 
the experimental results for the characteristic length. 
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III. FRACTIONAL REACTION-SUBDIFFUSION EQUATION WITH UNIFORM REACTIVITY 
A. Derivation from a mesoscopic CTRW model with reaction 

The starting point to derive the reaction-subdiffusion equation is the fundamental equation of a particle performing 
a continuous time random walk in the presence of a first order evanescence reaction (degradation). Consider the 
stochastic motion of an evanescent particle performing nearest-neighbor jumps on an infinite one-dimensional lattice. 
We shall hereafter denote the distribution of its waiting time between consecutive jumps by ^(^-The probability that 
a particle starting at site arrives at site i in the time interval between t and t + dt can be written as q(i, t) dt, where 
the arrival density q obeys the following integral difference equation: 



«(*.*) = f dt'ip RW (t-t')q(i + (-l)i,t') + S i)0 . (4) 

In this equation, ipRw(t) = ip(t)e~ kt is the probability per unit time that a particle found at a given site at time 
t = has performed a jump up to time t in the presence of a uniform evanescence reaction (since evanescence and 
jump are independent processes, the probability of jump is simply multiplied by the survival probability e~ k '). Taking 
the diffusive limit of Eq. Q one can show via suitable Fourier-Laplace techniques [f| that the probability w(x, t\0, 0) 
to find the particle at position x after a time t given that it was initially at x — obeys the following equation: 



=Kl e-- Vl-^^ W (x,t\0,0) 

-kw(x,t\0,0). (5) 
The operator o "D t 7 is defined via the equation 



C-}, t {u 1 -ry(u)}= o^yW, (6) 

where y{u) is the Laplace transform of the function y(t) and C~^ t {•} denotes the inverse Laplace transform. The 
operator o "D t ~ 7 is closely related to the Ricmann-Liouville fractional derivative 



D l t -y{ x ,t) = 



1 d 



dt' ■ 



f(x,t') 



r(i)dtj (t-t>) 



(7) 



In fact, qT>]~ j and oD\ ~ 7 are the same when applied to sufficiently regular functions f(t), as determined by the 
condition lim t _>.o J dt' ' (t — t') 7_1 /(£') = 0. This condition is actually fulfilled by all functions of t relevant to the 
morphogen problem, hence we shall use oD\ _1 in place of o^t" 7 in what follows. 

If one is dealing with more than one particle, the concentration c(x, t) follows the same kinetics as above, i.e. [1, Q 



dc(x, t) 
dt 



K 1 e 



-kt 



kc(x, t). 



(8) 



As one can see, this equation is a non-trivial extension of Eq. ([T]) for the case of anomalous subdiffusion with 
constant reactivity k. In the normal diffusion limit 7 — » 1 the Riemann-Liouville operator reduces to unit and one 
recovers Eq. ([T]) with a constant k. On the other hand, in the absence of reaction (k — > 0) Eq. ([5]) reduces to the 
standard fractional diffusion equation, which yields sublinear growth of (x 2 ). 

Turning now to the morphogen problem, Eq. ([8]) is to be solved subject to the boundary condition ([2]). The solution 
c{x 1 1) for the case of a particle source can be obtained from the propagator solution cp(x, t) = G(x, t) (corresponding 
the initial condition G(x,0) = S(x)) via the relation c(x, u) = jo G{x,u)/u between the Laplace transforms. The 
solution in Laplace space is found to be 



c{x,u) = — -= — exp 



-(u + kyi l /^JK^\x\ 



(9) 
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The stationary solution is obtained from the final value theorem for the Laplace transform: 

fc7 /a-i 



c s (x) = lim uc(x,u) 



Jo ■ 



u-yO 



IK, 



■ exp 



-\x\k^ /2 /^/KZ > 



(10) 



Eq. (fTUI) generalizes the exponential profile described by Eq. . Steady state profiles are thus seen to also exist 
in the presence of anomalous diffusion, as opposed to what had been suggested in some previous works [i[ . 



B. Robustness of stationary profiles 

Using Eq. (fTU)) it is possible to study the robustness of the concentration profiles with respect to a perturbation in 
the incoming flux jo. To this end, we take a reference value c x of the concentration and assess how large the shift of 
the associated position 



2c* ,/iL 



(11) 



becomes when jo is perturbed; the larger the shift, the smaller the robustness of the profile. The latter can thus be 
characterized by the inverse of the relative change of x with respect to a characteristic length a of the problem (e.g. 
the linear size of a cell), i.e. 



dx 



(12) 



Inserting Eq. (Ill[) into this definition we find 




(13) 



IV. FRACTIONAL REACTION-SUBDIFFUSION EQUATION WITH NON-UNIFORM REACTIVITY 

Seki et al. [l(| have shown that a CTRW process described by a generalization of Eq. ((4]) , namely 

1 rt 



3=0 
Si,0 



(14) 



with ipRw{i,t) = Tp{t)z yields the following reaction-subdiffusion equation: 



dc(x, t) 

at 



7 dx 2 

-k(x)c(x, t). 



(15) 



In order to tackle the corresponding morphogen problem, it is first necessary to find the propagator solution of Eq. 
T5|) . To this end it is convenient to introduce a new function v(x, t) defined via the transformation 



v(x, u) = [u + k(x)] 1 7 c(x, u) 
in Laplace space. This function is readily found to fulfil the equation 



(16) 



d 2 

[u + k{x)Vv(x, u) — 6(x) — Ky—^rv(x, u). 

ox 1 



(17) 



In what follows, Eq. (fTT|) will be used to investigate the effect of a non-uniform reactivity for several special cases. 
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A. Piecewise constant reactivity 

Here, we assume that the reactivity is given by a superposition of Heaviside functions, i.e. k(x) = koH(R — |x|) + 
kiH(\x\ — R). In region (0 < |x| < R) one has k(x) — k > 0, whereas in region 1 (|x| > R) one has k(x) = k\ > 0. 
Let us respectively denote by Vq(x, u) and Vi(x,u) the solutions of Eq. (IT71) in the regions and 1. These functions 
must fulfil the continuity conditions 



Vq(R,u) = Vi(R,u) 



and 



dva(x, u) 



dx 



dv\{x, u) 



\x\=R 



dx 



(18) 



(19) 



\x\=R 



In contrast, an integration of Eq. (|17[) across the origin shows that the solution must be discontinuous there: 



dvo(x, it) 



dx 



dvo(x, u) 



x=0+ 



Ox 



x=0- 



(20) 



Using Eqs. (I18[) - (I20D one can find explicit expressions for the Laplace transforms v(x,u), G(x,u) and c(x,u) 
j{)G{x,u)/u. For |x| < R one gets the stationary biexponential solution 



c s (x) = j kl 1 w (a;, u-i-0) 



(21) 



with 



V (x,u) = A e' aoX + B e ao:i 
The characteristic constants are 



(u + fco) 7 



(22) 



A 



Bn = 



{k + u)-^' 2 /2JK' 1 



2R(fc + T ,)T/ 2 

1 ) e +1 



-l-(k 1 +u)-i/ 2 (k +u)--i/ 

(ko + u)-y/ 2 /2^/K^ 

2R(k + u)~f/ 2 

--1 e -1 



l-(ko+u)--< / 2 (ki+u)~i 

For \x\ > R we shall distinguish two subcases with different physical behaviour. For k\ > one asymptotically gets 
the exponential decay law 



C s {x) = c(x,t -> GO) OC e -kl /2 (x-R)/^K 



(23) 



In contrast, for k\ = the behaviour is different. For normal diffusion the profile becomes constant for large |x| 



i.e., 



c s (x)<xj (h = 0,7 = 1, |x| > R). 
However, when the diffusion is anomalous one has 



(24) 



z(x, t oo) oc j t 1 7 (ki = 0, 7 < 1, |x| > i?), 



(25) 



i.e., there is no steady state! In view of Eq. (|2"Tj) and (|2"5j) . we conclude that the profile is discontinuous at x = R. 
This behaviour is confirmed by numerical simulations (see Fig. [1}. 
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FIG. 1: Simulation results (symbols) of c(x,t) for a step reactivity [k(x) = koH(R— x)] with fco = 1/1000 and R — 5 for 
7 = 1/2 (only values for x^O are shown). The particles are simulated by means of a CTRW model with the Pareto waiting 
time distribution l/>(t) = 7/(1 + t) 1+J and equiprobable jumps { — 1, 0, 1}. These parameters lead to the Kj-vahie \/\/Qk. The 
solid line corresponds to the theoretical prediction for the steady-state profile when x < R. For x > R no stationary profile is 
developed. The convergence of the simulation results to the stationary profile for x < R is very slow for values of x close to 
the discontinuity at x — R. No adjustable parameters were used. 

B. Exponentially decaying reactivity 

Here, we assume the decay law k{x) = ioe - ^' 1 '. While in this case Eq. (fl~7f does not seem exactly solvable for 
finite u, it is possible to find an exact expression of the steady state profile by techniques similar to the ones used 
above. The final result is 



ut/2-1 7, 



Cs{x) = jo- 



o (akl 



7/2 p -f3 1 \x\/2 



2K, 



1/2 



-(7-l)/8M 



(26) 



where the I„'s are modified Bessel functions and a = 2/ (/3"/WK^). As in the case of piecewise reactivity with 
k% = 0, this expression displays a different behaviour for normal and anomalous diffusion. In the normal diffusion 
case (7 = 1) one gets a monotonically decreasing profile from the concentration value 



c s (x = 0) 



3o 



(akl 



/2 



at the origin to the limiting value 



c s (x — ¥ ±00) 



Jo 



(27) 



h (akl /2 ) 



(28) 



(see Fig. [5]). In contrast, for 7 < 1 we find a qualitatively different behaviour. As one moves away from the source, 
first the concentration decreases until it reaches a minimum and then it increases (see Fig. [3]). 



C. Power law reactivity 



Next, we take k(x) = kq(xq + \x\) 1 with \x > 0. Since a steady state was already attained for an exponen- 
tially decaying reactivity, this will also be the case under the present situation, which describes enhanced particle 
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FIG. 2: Convergence of CTRW simulation results (symbols) to the stationary profile predicted by formula (|26[l for jo = 1,7 = 1, 
the corresponding value of the diffusion coefficient K 1 = 1/3 and exponentially decaying reactivity k(x) = fco exp(— f3\x\) with 
fco = 1/100 and /3 = 1/8 (solid line). CTRW jump characteristics as in fig. [1] 




FIG. 3: Convergence of CTRW simulation results (symbols) to the stationary profile predicted by the formula (I26[) for jo = 1, 
7 = 0.5 , the corresponding value K-y = 1/V97T and exponentially decaying reactivity k(x) = fco exp(— ft\x\) with fco = 1/200 
and f3 = 0.6 (solid line). CTRW jump characteristics as in fig. [TJ The simulation results clearly go towards the stationary 
solution as time increases, although the convergence for large x is slow. 



evanescence. The general solution of Eq. (TP7|) for fi ^ 2 is given by the modified Bessel functions Ji„i and Ki u i with 
v = (/i — 2) _1 . In order to single out the Bessel function corresponding to the physical solution, we use the fact that 
the incoming flux must be equal to the amount of particles per unit time that disappear due to degradation, i.e., 

/oo 
k(x)c s (x)dx. (29) 
-oo 

This condition leads to different solutions depending on the value of /j,. For fj, < 2 one gets 

c s (x) cx (x + \x\)^-^K M ($(x + M)-&) , (30) 

where $ = 2\v\^ k^/ K^. For large \x\, the above stationary solution can be shown to go to zero as 
x^~~ exp f— <&\x\ In contrast, for /i > 2 one has 
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c s (x) oc (x + |a;|)*- /i+ '7 M ($(x + M)~*) • (31) 

As |x| — > oo, this expression tends to a constant limiting value in the normal diffusion case and grows as x~~^ for 
7 < 1. 

When [i = 2 the solution is not given by a power law rather than by Bessel functions. One has 

c s (x) oc (x Q +x)^ 2+X - (32) 

with A_ = 1 ± yl + 4(kq/A' 7 )/2. For |a;| — > oo this solution goes to infinity, a constant value or zero depending 
on whether — — 2 + A_ is positive, zero or negative. 

V. CONCLUSIONS AND OUTLOOK 

In the present work we investigate both analytically and numerically the behaviour of the stationary concentration 
profiles arising from fractional reaction-subdiffusion equations derived from a CTRW model with a superimposed 
death process. These fractional equations are a natural extension of classical reaction-diffusion equations traditionally 
used to study the problem of morphogen gradient formation. We consider the case of linear degradation with both a 
uniform and non-uniform reactivity. The formulation of the problem in terms of fractional diffusion equations turns 
out to be a key ingredient in the analysis of the properties of morphogen gradients. This approach allows us to exploit 
a plethora of powerful analytical techniques available from fractional calculus to tackle the morphogen problem. 

In the uniform case one obtains exponentially decaying stationary concentration profiles. Their robustness with 
respect to changes in the incoming flux increases with increasing k. Likewise, one can study the robustness of the 
profiles with respect to a perturbation in k by introducing the quantity 



U k = a (k-) (33) 
This issue will be the subject of future research. 

In the non-uniform case the behaviour of the stationary profiles turns out to be very sensitive to the specific spatial 
dependence prescribed for k and to the value of the anomalous diffusion coefficient 7. Moreover, for the case of a 
piecewise constant reactivity with k\ = and anomalous diffusion, we see that a discontinuous profile arises and no 
steady state is reached in region 1. This is a novel effect not seen for normal diffusion. For exponentially decaying 
reactivity the concentration goes to a constant limiting value far away from the source (\x\ — > 00) when 7 = 1, but it 
grows without bound for 7 < 1. Finally, when the connectivity decays as a power law, the stationary concentration 
may go to zero, to a constant or to infinity depending on the values of the characteristic decay exponent and 7. 

In view of the strong inhomogeneities encountered by the diffusing morphogens in the embryonic environment, we 
believe that the sensitivity of the concentration gradients to the form of k(x) may be relevant for the modeling of 
morphogen gradient formation and interpretation. 

Up to the case of piecewise constant reactivity with ki = 0, in the present work we limit ourselves to study the 
behaviour of the stationary concentration profiles. However, analytic solutions for transient profiles are available 
for some of the cases studied, and others can be investigated via numerical techniques for the inversion of Laplace 
transforms. Besides, a reaction-subdiffusion equation was recently obtained that generalizes the above results to the 
general case k = k{x,t) [O, HI]- As one could have guessed in view of Eqs. (JSj) and (fT5l) . this equation reads as 



dc(x, t) 
dt 



dx- 



{(>-/«? oiji-r Ut^'W' c[x,t)] \ - k(x,t)c(x,t) 



(34) 



Our aim is to use the above equation to investigate further problems related to morphogen gradient formation in 
future. Beyond this field of research, Eq. (f34|) can be applied to many other problems of interest characterized by 
different kinds of boundary conditions. 
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